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Abstract 


Overpressure  time  history  data  from  warhead  blast  experiments  yield  peak  overpressure  P 
as  a  function  of  spatial  position. 

Dr.  Owen  Litt  has  proposed  a  model  for  P  based  on  the  peak-overpressure  characteristics  of 
a  bare  spherical  charge.  The  direction-independent  peak-overpressure  function  for  a  bare 
spherical  charge  is  modified  to  have  nonspherical  level-surface  structure  by  specifying  surfaces 
of  constant  peak  overpressure.  This  introduces  a  directional  component  into  the  model. 

In  this  report,  the  original  formulation  is  refined  and  generalized  and  a  mathematical  model 
and  computer  code  are  presented  to  evaluate  the  function.  Such  a  computational  device  is 
required  for  model  parameter  estimation  and  experiment  design. 
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1.  Background 


Overpressure  time  history  data  from  warhead  blast  experiments  yield  peak 
overpressure  P  as  a  function  of  spatial  position.  Dr.  Owen  Litt  [1]  has  proposed  a  model 
for  P  based  on  the  peak-overpressure  characteristics  of  a  bare  spherical  charge.  In  that 
model,  the  direction-independent  peak-overpressure  function  for  a  bare  spherical  charge 
is  modified  to  have  a  level  curve  structure  with  a  specific  nonspherical  functional  form. 
This  induces  a  directional  component  into  the  model. 

It  is  necessary  to  combine  the  peak-overpressure  function  representations  of  the 
bare  spherical  charge  and  an  arbitrary  level-curve  structure  to  produce  the  required 
mathematical  model.  This  report  details  the  solution  of  that  analytical  problem  and  also 
an  explicit  solution  to  the  problem  of  level-curve  model  specification  in  general,  and  so 
serves  a  twofold  purpose.  The  development  of  a  general  theoretical  framework  for  solving 
such  model-specification  problems  appears  in  section  2.  The  rest  of  this  report  describes 
the  application  of  the  general  principle  to  the  specific  problem  of  Dr.  Litt’s  blast  model.  The 
original  formulation  is  refined  and  generalized  and  a  mathematical  model  and  computer 
code  are  presented  to  evaluate  the  function. 


2.  Level-Curve  Model  Specification 


This  section  contains  a  discussion  of  model  formulation  based  on  specifying  a  geometric 
level-curve  structure.  An  example  precedes  the  development  of  a  general  method  for  the 
formulation  of  such  models. 


2.1  An  Example.  Suppose  F  is  a  decreasing  function  on  [0,  °°).  For  example,  take 
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W=l+rr 

The  function  F  can  be  used  to  create  function  Fi :  R2  — ♦  R+  by  defining 


(1) 


FiM)  =  F(r)  =  — (2) 

where  the  usual  Cartesian  coordinates  on  K2  are  (x,y)  and  polar  coordinates  (r,</>)  on 
R2  are  given  by  x  =  rcos<j>  and  y  =  rsin<j>.  In  the  x-y  plane,  level  curves  of  Fi  are 
concentric  circles  centered  at  the  origin,  since  Fx  is  independent  of  <j>.  For  any  L  €  (0, 1], 
the  value  of  r  that  makes  Fi(r,  <f>)  =  L  is  given  by  Fi(r,  0)  =  F (r)  =  1/(1  +  r2)  =  L,  so  that 
r  =  F ~ 1  (I )  =  y/l/L-1.  In  other  words,  F ( y'l/L-l)  =  L ,  so  Fi  =  I  on  the  circumference 
of  a  circle  with  radius  i/l/L-1  and  area  n  ( 1  /L  - 1 ) .  It  is  possible  to  modify  the  definition 
of  Fi  and  produce  a  function  F2  that  has  elliptical  level  curves  with  a  given  eccentricity  and 
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orientation.  Furthermore,  the  value  of  Fi  on  a  given  circle  should  be  equal  to  the  value  of 
F2  on  an  ellipse  with  the  same  area.  The  polar  equation  for  an  ellipse  is 


r2  =  ab  •  \/l  -  e2  • 


1+tan  2(0-(/)o) 

1  -£2  +  tan2(0-0o)’ 


(3) 


where  the  ellipse  has  eccentricity  e  =  y/l  -b2/a2,  major  axis  length  2 a  in  the  direction 
<p  =  <t>0,  minor  axis  length  2b  in  the  direction  </>  =  4>0  +  n/2,  and  area  nab.  The  function  F2 
is  now  defined  by 


F2(r,</))  =  F 


l  +  tan2(<ft  — 0O) 

1  -e2  +  tan2(0-0o)J 


1-1/2' 


(4) 


Then,  it  can  be  seen  that  F2(r,  <t>)=L  when 


r 
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1  +tan2(<ft-<ft0) 

1  -£2+tan2(0-0o)’ 


(5) 


so  that  F2(r,  <p)  =  L  on  the  perimeter  of  an  ellipse  of  area  tt(1  /L  - 1).  (See  Figures  1  and 
2  for  a  depiction  of  the  level  curves  of  these  example  functions.) 


2.2  The  General  Construction.  Consider  a  decreasing  function  F:  [0,°°)  -*  R+. 
This  function  F  can  be  used  to  create  a  function  Fi :  R2  — »•  R+  by  defining 

F1(r,0)  =  F(r),  (6) 

where  r  and  <f>  are  polar  coordinates.  Since  Fi  is  independent  of  <p,  the  level  curves  of  Fi 
are  concentric  circles  centered  at  the  origin  of  R2.  And  because  F  is  decreasing,  the  value 
of  Fi  is  smaller  on  a  larger  such  circle. 

It  is  possible  to  construct  a  version  of  F  that  has  noncircular  level  curves  with  any 
specific  functional  form.  In  particular,  say  the  level  curves  are  to  be  given  by 

r(4)=g*(u)  (7) 

for  various  values  of  u.  Suppose  for  all  <j>  that  g^(u)  is  a  continuous  function  of  u,  that 
g*(0)  =  0,  that  g*(u)  is  an  increasing  function  of  u,  that  g<p(u)  is  defined  for  all  u  ^  0,  and 
that  limu_oog^(u)  =  So,  g0  is  a  bijective  function  on  [0,°°),  and  g^  has  an  inverse  in 
the  following  sense:  for  each  fixed  value  of  <j>,  the  inverse  function  g^1,  characterized  by 
g«  (g*->))  =g$1  (g<p(u))  =  u,  is  well-defined  for  u  ^  0.  In  fact,  g^(u)  is  also  an  increasing 
function  of  u  on  [0,  oo). 

A  function  F2 :  R2  — ►  R+  satisfying  equation  (7)  can  be  defined  in  terms  of  polar 
coordinates  by 

h  (r,<P)  =  F(g;'(r)).  (8) 

To  show  this,  let  u  >  0  be  constant.  Then  F  (u)  is  also  constant,  and  the  locus  of  (r,  <p)  which 
has  F2(r,0)  =  F(u)  is  given  by  F(u)  =  F2(r,<p)  =  F  (g^(r)).  This  means  that  u  =g"1(r), 
from  which  equation  (7)  follows. 
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3.  Application  to  Blast  Model  Formulation 

Here,  the  results  of  the  previous  section  are  applied  to  the  formulation  of  models  for  the 
maximum  peak  overpressure  of  a  detonating  charge  blast  field.  First,  the  basic  formulation 
is  discussed  and  then  an  enhanced  model  is  presented. 


3.1  Basic  Model.  This  model  works  in  two-dimensional  polar  coordinates  (r,  <f>)  with 
the  origin  centered  on  the  detonating  charge.  A  three-dimensional  spatial  model  for  peak 
overpressure  P  can  be  obtained  by  rotating  a  two-dimensional  model,  defined  in  the 
half-plane  0  <p  <:  n,  about  the  x-axis. 

The  development  of  a  two-dimensional  model  for  peak  overpressure  as  a  function  of 
the  polar  coordinates  (r,  <j>)  follows.  In  this  report,  models  for  peak  overpressure  are  based 
on  the  function  Ps(z),  which  gives  the  maximum  peak  overpressure  for  detonation  of  a 
spherical  TNT  charge.  In  the  definition  of  Ps  and  throughout  this  report,  the  normalized 
distance  coordinate 


is  used,  where  r  is  measured  in  feet,  W  is  charge  weight  in  pounds,  and  a  is  a  constant 
with  nominal  value  a  =  1/3.  The  spherical  charge  pressure  function  Ps  ( z )  itself  is  defined 
by 


Ps(z)  =  exp 


A 

z  +  B 


(10) 


where  the  constants  have  the  empirical  values  A  =  31.97,  B  =  3.555,  and  C  =  0.5.  This 
function  was  derived  from  a  fit  to  empirical  data  [1]. 


The  modeling  concept  under  consideration  requires  a  pressure  function  component  Pn 
with  a  nontrivial  dependence  on  <p,  specified  by  a  certain  family  of  noncircular  level  curves. 
The  function  P„  is  derived  from  Ps  in  the  same  way  that  F2  is  derived  from  F  in  section  2,  by 
the  application  of  equation  (8)  to  a  specific  level  curve  function.  The  level  curve  function 
for  Pn,  which  gives  an  appropriate  shape  based  on  engineering  considerations,  is  given  by 


g*(u)  =  u 


4  n(u)/m 

5 


where  n(u)  is  defined  shortly,  and  Pn  is  defined  by 

Pn(z,<t>)=Ps{g;1(z)), 


(ID 


(12) 


as  in  section  2. 
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Two  more  definitions  complete  the  specification  of  Pn.  The  function  /  is  defined  to 
be  a  “smooth  step  function”  with  /( 0)  ~  0,  /  increasing,  and  /(z)  — ►  1  as  z  — ►  °°.  To  be 
specific,  set  z0  =  10  and  take 

/(z)  =  i  +  iarctan(z-20).  (13) 

z  7T 

The  function  n  is  defined  by 


n(u)  =  n0(l-/(u)),  (14) 

where  n0  is  a  positive  constant.  So  n  is  a  decreasing  function  with  n(0)  =  n0  and  n(u)  — ►  0 
as  u  — ►  oo.  The  behavior  of  n  along  with  the  form  of  g#  implement  the  design  objective 
that  Pn  looks  like  Ps  at  large  distances;  i.e.,  the  level  curves  of  Pn  become  circular  for  large 
z. 


Now  with  the  function  Pn  completely  specified,  the  conditions  on  g^(u)  of  section  2 
are  indeed  satisfied.  The  exponent  4 n(u)/m  is  positive  and  decreasing  with  u,  so 
sin(m0/2)4n(u)/m  is  a  nondecreasing  function  of  u  for  fixed  <f>.  Therefore,  g^(u)  is 
increasing  in  u  for  any  0.  The  conditions  of  section  2  are  satisfied,  so  a  level  curve  of 
Pn(z,(p)  is  given  by  z  =g^(u)  for  u  fixed,  as  required. 

To  specify  the  peak  overpressure  model,  it  remains  only  to  combine  the  function 
components  Ps  and  Pn  in  a  certain  way.  The  definition  of  the  peak  overpressure  model 
function  P  is 


P(z,0)  =  /(z)Ps(z)  +  (l-/(z))Pn(z,0),  (15) 

where  the  functions  Ps,  Pn,  and  /  are  as  previously  discussed.  Because  of  the  nature  of 
/,  the  pressure  function  looks  like  Pn  for  small  z  and  like  Ps  for  large  z.  Definition  of  the 
model  is  now  complete.  The  quantities  P,  z,  /,  P$,  g#,  and  n  were  specified  by  Dr.  Owen 
Litt  [1,  2,  3,  4,  5],  as  was  an  implicit  characterization  of  Pn.  The  explicit  representation  of 
equation  (12)  for  Pn  is  a  product  of  this  report. 

In  summary,  the  complete  model  is  given  by 

P(z,<p)  =/(z)Ps(z)  +  (l-/(z))Pn(z,</>),  where 
z  =  r/Wa, 

f(z)  =  l/2  +  l/7r-arctan(z-z0), 

P*(z)  =  exp(A/(z+B)-C), 
n(u)  =  n0(l-/(u)), 
g*(u)  =  u  (sin  m<f>/2)4n{u)/m ,  and 
Pn(z,<p)  =  Ps(g<p-1(z)).  (16) 

The  quantities  A,  B,  C,  and  W  are  constant;  z0  =  10,  Zi  =  5,  and  z2  =  15  are  fixed  model 
parameters;  and  m,  n0,  and  a  are  model  parameters  to  be  estimated.  Interpretations  of 
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the  parameters  are  as  follows:  m  determines  the  direction  of  the  Pn  component,  the  value 
of  n0  makes  the  Pn  component  more  or  less  concentrated  in  the  direction  determined  by 
m,  and  a  determines  the  dependence  of  normalized  distance  on  charge  weight. 

Now  some  characteristics  of  the  model  can  be  examined  in  more  detail.  The  extreme 
point  on  a  level  curve  occurs  when  <p  =  n/m  and  sin(m0/2)  =  1,  in  which  case  gn/m(u)  —  u 
and  also  g~}m(u)  =  u.  Then,  Pn(z,n/m)  =  Ps(z).  So,  in  the  direction  of  maximum  peak 
overpressure,  <p  =  n/m,  the  Pn  component  has  the  same  pressure  value  as  spherical  bare 
charge,  Ps.  In  other  directions,  for  fixed  z,  the  value  of  Pn  is  lower  than  Ps. 

For  an  example,  set  the  charge  weight  to  W  =  1  and  set  the  function  parameters  to 
m  =  1.75,  n0  =  2.0,  and  a  =  1/3.  Figure  3  demonstrates  the  level  curve  characterizations 
of  Pn  and  Ps  at  the  same  pressure  value.  Pn  is  evaluated  at  the  point  (z0,  <p0)  =  (4,  n/3). 
This  point  lies  on  the  level  curve  z  =g(t>(u),  where  g^  (u)  =  z0,  or  u  =g^0~1(z0),  so  a  general 
point  on  this  level  curve  has  coordinates  (g<p(u),  <p).  The  extreme  point  on  this  level  curve, 
where  <p  =  n/m,  has  coordinates  (gn/m(u),  n/m)  =  (u,  n/m).  The  value  of  Pn  anywhere 
on  this  level  curve  is  Pn((g*(u),  0))  =  Ps(g^\g^)))  =  Ps(u).  Particular  values  for  this 
example  are  u  ~  8.94  and  Ps(u)  ^  8.49. 

Once  again,  in  the  direction  of  maximum  peak  overpressure,  0  =  n/m,  the  P„ 
component  has  the  same  pressure  value  as  spherical  bare  charge,  Ps.  In  other  directions, 
for  fixed  z,  the  value  of  Pn  is  lower  than  Ps .  This  may  not  be  realistic.  In  the  maximum 
direction,  Pn  should  have  a  higher  value  than  Ps,  since  the  blast  modeled  by  P„  is  focused  in 
that  direction.  This  additional  feature  is  implemented  by  incorporating  into  the  definition 
of  Pn  an  equivalent  spherical  charge  weight  Ws  that  is  greater  than  the  actual  charge 
weight  W,  effectively  renormalizing  the  distance  coordinate  in  equation  (9),  which  is 
then  used  in  equations  (10)  and  (12).  A  conceptually  equivalent  approach  is  to  directly 
reduce  the  distance  argument  z  of  Ps  in  equation  (10),  as  it  is  used  in  the  definition  of 
P„.  Alternatively,  the  level  curve  function  can  be  changed  in  the  definition  of  P„  from  g^, 
equation  (11),  to  a  new  function  produces  a  higher  pressure  value  on  the  level  curve. 

As  shown  in  section  3.2,  these  three  schemes  are  equivalent.  The  net  effect  of  any  of 
them  is  to  force  a  P„  level  curve  to  correspond  to  a  smaller  Ps  level  curve,  on  which  the 
pressure  is  higher.  The  basic  model  can  be  modified  to  have  this  property. 


3.2  Enhanced  Model.  The  model  of  the  section  3.1  is  generalized  by  the 
introduction  of  a  new  pressure  function  P*,  which  replaces  Pn.  The  P*  component  in 
the  maximum  direction  (0  =  n/m)  has  the  pressure  value  of  a  spherical  charge  of  arbitrary 
weight  Ws(z).  To  increase  the  generality  and  flexibility  of  the  model,  Ws  is  allowed  to  be 
a  function  of  z  rather  than  a  constant.  It  is  convenient  to  define  the  function  M  by 


M(z)  = 


Ws(z) 
W  ’ 


(17) 


so  that  M  represents  a  dimensionless  mass  scaling  ratio  or  magnification  factor  in  the 
maximum  direction,  since  Ws(z)  is  the  equivalent  sphere  charge  weight  in  that  direction 
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(18) 


(n/m)  at  the  distance  r  =  zWa.  Now,  define  P*  by 

Pn*(z,0)  =  Ps(s(g^(z))) 

and  proceed  to  solve  for  s.  When  $  =  n/m,  the  result  is 

P*(z,n/m)  =  Ps  (s(g-}Jz)j)  =  Ps{s(z)). 

Since  z  =  rW~a  and  s(z )  =  rWs(z)_a  =  rW~aM(z)~a, 
r  =  zWa  =s(z)WaM(z)a,  and  then 


M(z)  = 


or  s(z)  =zM(z) 


(19) 

it  follows  that 


(20) 


This  expresses  the  required  function  s  in  terms  of  the  magnification  factor  M  and, 
therefore,  also  in  terms  of  the  equivalent  sphere  weight  Ws.  Note,  referring  to 
equation  (10),  that  the  function  s  as  it  appears  in  Ps(s (z))  amounts  to  a  rescaling 
of  distance  in  the  function  Ps.  Also,  the  definition  of  P*  can  be  written  as 
P*(z,<j>)  =  Ps((g<p^h)~1(z)),  where  ft-1  =s,  and  thus  making  explicit  the  modification  of  the 
level  curve  function  in  the  definition  of  P*.  So  the  three  conceptual  approaches  (weight 
scaling,  distance  scaling,  and  level  curve  modification)  to  the  derivation  of  P*  from  Pn  are 
equivalent. 

Now,  with  Ws(z)  =  W,  then  M(z)  =  1  and  s(u)  =  u,  so  that  P*  =  Pn.  This  reduces  to 
the  basic  model  of  the  section  3.1,  where  the  Pn  has  the  property  that,  in  the  maximum 
pressure  direction  <f>  =  n/m,  the  peak  overpressure  is  equal  to  that  of  a  spherical  charge  of 
the  same  weight. 

A  more  realistic  general  formulation  requires  that  M(0)  >  1,  that  M(z)  is  a 
nonincreasing  function  of  z,  and  that  M(z)  — ►  1  as  z  — ►  °°.  This  behavior  embodies  the 
design  criteria  that  P*  itself  looks  like  Ps  at  large  distances  and  that  the  pressure  P*  is 
greater  than  Ps  in  the  maximum  direction  at  small  distances. 

It  may  be  possible  to  completely  specify  M  through  energy  conservation  considerations, 
but,  for  illustrative  purposes,  a  piecewise  continuous  version  of  M  is  used.  This  has 
corresponding  s  that  is  easy  to  calculate.  Let  M(z)  =  M0  >  1  for  z  <  zj,  let  M(z)  =  1 
for  z  >  z2,  and  let  M(z)a  be  a  linear  function  of  z  for  zj  <  z  <  z2.  It  is  convenient  to 
express  piecewise  function  definitions  in  terms  of  the  indicator  function 


IT(t)  = 


t  e  T 
t&T. 


(21) 


First  define  the  “linear  step  function”  L  with  the  characteristics  that  L(z)  =  b  for  z  ^  z1; 
L(z)  =  1  for  z  ^  z2,  L(z )  is  linear  for  Zi  <  z  ^  z2,  and  L  is  continuous.  The  appropriate 
definition  is 


L(z;  b, zi, z2)  =  b  ■  J[0,Zl)(z)  +  (ai  +a2z)  •  J[Zl,22](z)  + 1  •  J(Z2,„)(z),  (22) 
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bz2-zx 

where  a.\  = - 

z2-zi 


1-5 


(23) 


and  a2  = - . 

z2-zi 


The  corresponding  definition  for  M  is  then 

M(z)=L(z;M;,z1,z2)1'‘.  (24) 

Solve  for  5  in  closed  form  to  get 


s(z)  =  z/b  ■  I[0iZl)(z)  +z/(d!  +a2z)  ■  IlzuZ2](z)  +z  ■  I(Z2iCc}(z),  (25) 

where  b  =  M“,  and  ai  and  a2  are  given  by  equation  (23). 

A  complete  working  model  is  then 

P(z,<j>)=f(z)Ps(z)  +  (l-f(z))p*(z,<p)  where 
z  =  r/Wa, 

f(z )  =  l/2  +  l/7r-arctan(z-z0), 

Ps(z)  =  exp (A/(z  +  B)  -  C), 
n(u)  =  n0(l-/(u)), 
g*(u)  =  u  (sin  m<p/ 2)4n(u)/m , 

M(z)  =  L(z;  M“,  Zl,  z2)ly/a, 
s(z)  =  zM  (z)~a,  and 

p;(z,4,)  =  ps(s(g;'(z))).  (26) 

The  quantities  A,  B,  C,  and  W  are  constant;  z0  =  10,  z  1  =  5,  and  z2  =  15  are  fixed  model 
parameters;  and  m,  n0)  M0,  and  a  are  model  parameters  to  be  estimated. 

The  example  of  section  3.1  illustrates  the  enhanced  model.  Again,  the  charge  weight 
is  W  =  1  and  set  the  function  parameters  are  m  =  1.75,  n0  —  2.0,  and  a  =  1/3.  The  new 
function  parameter  for  mass  scaling  is  M0  =  4.  Figure  4  depicts  the  functions  M  and  s. 
Figure  5  demonstrates  the  level  curve  characterization  of  P*  in  relation  to  that  of  Ps  at  the 
same  pressure.  As  before,  P*  is  evaluated  at  the  point  (z0,  <p0)  =  (4,  n/3).  This  point  lies 
on  the  level  curve  z  =g</,(u)  where  g^0(u)  =z0,  or  u  =g<p0~1(z0),  so  a  general  point  on  this 
level  curve  has  locus  (g^(u),  </>).  The  extreme  point  on  this  level  curve,  where  <p  =  rc/m, 
has  coordinates  (g„/m(u),  n/m)  =  (u,n/m). 

Note  that  the  P*  level  curve  is  identical  to  the  Pn  level  curve  in  Figure  3,  which 
illustrates  the  example  of  section  3.1.  The  function  value  is  different,  however, 
to  reflect  the  increased  equivalent  sphere  charge  weight  or  reduced  distance  in  Ps. 
The  corresponding  Ps  level  curve  in  Figure  5  has  a  radius  smaller  than  the  extreme 
distance  on  the  P*  level  curve.  The  value  of  P*  anywhere  on  its  level  curve  is 
pz((g*(u) ’<!>))  =ps(s(gt1(g*(u))))  =  ps(s{u)),  which  is  also  the  value  of  Ps  on  its  level 
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curve  in  Figure  5.  Particular  values  for  this  example  are  u  ~  8.94,  s(u)  ~  6.59,  and 
Ps(s(u))  ~  15.6.  The  weight  scaling  factor  is 


M(u)  = 


u 


1/a 


~  2.49, 


(27) 


and,  since  W  =  1,  this  is  also  the  equivalent  sphere  charge  weight  in  the  direction  </>  =  n/m 
at  the  distance  r  =  u.  These  values  of  u,  s(u),  and  M (u)  are  distinguished  in  Figure  4. 


4.  Model  Evaluation  for  Experiment  Design 

Conducting  an  experiment  to  calibrate  the  model  (estimate  the  parameters)  involves 
placing  pressure  sensors  in  the  detonation  field  of  an  explosive  charge.  The  sensors  must 
be  placed  so  that  optimal  useful  information  is  obtained  from  the  experiment.  Sensors 
cannot  be  overdriven.  On  the  other  hand,  each  sensor  has  a  lower  limit  of  resolution, 
beyond  which  the  noise  in  the  measurement  system  overrides  any  signal.  Sensors  must 
also  be  placed  so  that  they  register  the  nonspherical  Pn  component  of  the  pressure  field. 
It  is  therefore  reasonable  to  “guess”  what  the  model  parameters  are,  evaluate  P,  and  place 
the  sensors  accordingly. 

A  graphical  display  of  the  model  response  is  useful  in  the  design  of  an  experiment  for 
blast  model  parameter  estimation.  Since  the  model  is  a  well-defined  function,  evaluation 
is  conceptually  simple:  replace  constants  with  numbers  and  evaluate  the  functions. 

The  form  chosen  for  M(z)  yields  a  closed-form  representation  for  s.  Generally,  g^1 
must  be  evaluated  numerically,  even  if  s  has  a  closed-form  representation.  Choosing 
another  form  for  M  ( z )  may  result  in  s  having  no  closed-form  representation,  which  will 
increase  computational  complexity. 

Figures  6-10  are  contour  representations  of  P  computed  with  parameter  values  W  =  1, 
m  =  1.75,  n0  —  2.0,  M0  =  4,  and  a  =  1/3  on  an  various  x-y  grids.  Spatial  coordinates 
are  equivalent  to  z  units,  since  W  =  1.  Contour  levels  are  indicated  in  the  figure  captions. 
Note  that  the  logarithmic  spacing  of  the  level  curves  gives  a  better  visual  display  than 
linear  spacing  would.  Figure  6  has  -20  x  <  20  and  0  ^  y  ^  20  to  show  the  far  field, 
Figure  10  has  —  1  ^  x  <  1  and  0  <  y  <  1  to  show  the  near  field,  and  the  intervening  figures 
depict  intermediate  ranges.  Computations  and  graphics  were  done  with  Mathematica  [6] ; 
the  code  necessary  to  reproduce  these  calculations  are  presented  in  the  Appendix. 


5.  Model  Parameter  Estimation 

Data  consist  of  empirical  measurements  of  peak  overpressure  p  at  spatial  location  (r,  0), 
denoted  as  p*,  r;,  and  0,  for  1  <  i  <  N.  Parameter  estimates  can  be  obtained,  for  example, 
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by  least  squares  in  the  response  or  log  response;  i.e., 

N 

minimize^  [p*  -  P(W~a ru  fa)]2,  (28) 

t=i 

or 

N 

minimize y^[logpi  -logP(W~arj,  fa)]2 ,  (29) 

i=  1 

where  the  minimizations  are  conducted  over  the  parameter  vector  (m,  n0,  M0,  a).  Due  to 
the  exponential  nature  of  P  and  the  error  characteristics  of  pressure  sensors,  estimation 
based  on  logP  will  most  likely  yield  more  accurate  results  than  estimation  based  on  P. 
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Figure  1.  Level  Curves  of  Fi  From  Section  1. 


Figure  2.  Level  Curves  of  F2  From  Section  1. 
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Figure  3.  Level  Curves  for  Pressure  Components  Ps  and  P, 


M(z ) 


s(z) 


Figure  4.  Mass  Scaling  Functions  M  and  s. 
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Appendix:  Mathematica  Code 
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This  following  Mathematica  code  is  provided  for  function  evaluation  and  visualization. 
This  environment  is  useful  for  preliminary  investigation,  function  selection,  and 
experiment  design.  Estimation  procedures  are  not  provided.  Mathematica  names  are 
generally  consistent  with  names  in  the  body  of  this  report. 


A.1  Evaluation.  Define  the  utility  functions  Seq  and  Lseq  to  create  linear  and 
logarithmic  sequences.  The  resulting  sequences  range  from  a  to  b  and  contain  n  elements. 

(*  UTILITY  FUNCTIONS  *) 

Clear [Seq,  Lseq] ; 

Seq [ a_ ,  b_,  n_]  :=  Range[a,  b,  (b-a)/(n-l) ] ; 

Lseq[a_,  b_,  n_]  :=  Exp [Seq [Log [a] ,  Log[b],  n] ] ; 

Define  the  model  constants  and  parameters  Ao,  Bo,  Co,  w,  alpha,  m,  No,  Mo,  zO,  zl,  and 
z2. 


(*  PARAMETERS  &  CONSTANTS  *) 

Clear[Ao,  Bo,  Co,  w,  alpha,  m,  No,  Mo,  zO,  zl,  z2]; 
(*  sphere  function  parameters  *) 

Ao  =  32.97; 

Bo  =  3.555; 

Co  =  0.5; 

(*  charge  weight  and  scaling  exponent  *) 
w  =  1.0; 
alpha  =  1/3 . 0 ; 

(*  model  parameters  *) 
m  =  7/4; 

No  =  2.0; 

Mo  =  4  ; 

(*  transition  function  parameters  *) 
zO  =  10;  zl  =  5;  z2  =  15; 

Define  the  model  functions  f  01,  f,  n,  k,  g,  gi,  Psz,  and  P. 

(*  BLAST  FUNCTIONS  *) 

Clear[f01,  n,  f,  s,  g,  gi,  Psz,  P] ; 

(*  basic  transition  function  *) 
f01[x_]  :=  1/2  +  ArcTan[x]/Pi; 

(*  generic  transition  fucntion  *) 
f[z_,  m_,  r_]  :=  f01[r(z-m)]; 

(*  level  curve  exponent  function  *) 
n[z_,  No_]  ;=  No  (  1  -  f[z,  zO,  1]); 


23 


(*  equivalent  sphere  weight  function  *) 
s[z_,  b_,  zl_,  z2_]  :=  z/b  /;  z  <=  zl; 

s [ z_,  b_,  zl_,  z2_]  := 

z/( (b  z2-zl)/(z2-zl)  +  (l-b)/(z2-zl)z)  /;  zl<z  &&  z<z2; 
s [z _ ,  b _ ,  zl _ ,  z2 _ ]  : -  z  / ;  z  >—  z2; 

(*  level  curve  function  *) 

g[z_,  phi_,  m_,  No_,  Mo_]  :=  z  Sin[m  phi/2] "(4  n[z,No]/m); 

(*  level  curve  inverse  function  *) 

gi[z_,  phi_,  m_,  No_,  Mo_]  :=  Module[{uO,  ul,  u}, 

If[g[z,  phi,  m,  No,  Mo]>z, 

For[uO=z,  g[uO,phi,m,No,Mo]  >  z,  u0=u0/3] ;  ul=3  uO , 
For[ul=z,  g[ul,phi,m,No,Mo]  <  z,  ul=3  ul] ;  u0=ul/3]; 
u  =  FindRoot [g[u,  phi,  m,  No,  Mo]  =  =  z,  {u,  {uO,  ul]], 
Maxlterations  ->  25]  [ [1] ] [ [2]  ]  ; 
u]  ; 


(*  sphere  charge  function  *) 

Psz[z_]  :=  Exp[Ao  /  (  z  +  Bo  )  -  Co]; 


(*  peak  overpressure  model  function  *) 

P[x_,  y_]  :=  Module [{r,  phi,  z,  uO ,  ul,  u,  P], 

r  =  Sqrt [x~2+y~2] ; 
phi  =  ArcTan[x,  y] ; 
z  =  r/w~alpha; 
u  =  gi[z,  phi,  m,  No,  Mo]; 

P  =  f [z,  zO,  1] Psz [z]  +  (l-f [z, zO , 1] ) Psz [s [u,Mo~alpha, zl, z2] ] ; 
P] 


A.2  Visualization.  Define  the  function  Ptable  to  evaluate  P (x,  y)  on  an  nx  by  ny 
grid  with  -XI  ^  x  <  XI  and  dy  ^  y  <  Xl+dy. 


Ptable [Xl_,  nn_]  :=  Module [ 

(nx,  ny,  dy  =  0.05,  X0,  Y0,  X,  Y,  Yl], 

X0  =  -XI;  Y0  =  0;  Yl  =  Xl ;  nx  =  nn;  ny  =  nn/2 ; 

X  =  Seq [X0 ,  XI,  nx] ; 

Y  =  dy+Seq [ Y0 ,  Yl,  ny] ; 

XP  =  Table[P[X[ [i] ] ,  Y[[j]]],  {i,  nx] ,  {j,  ny]]; 

{X0,  XI,  nX,  Y0,  Yl,  ny,  XP}] 
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Define  the  function  Pshow  to  graph  xp,  the  result  of  Ptable.  The  other  arguments  are 
the  lowest  contour  level  (LO),  the  highest  contour  level  (LI),  and  the  number  of  countour 
levels  (nl). 


Pshow [XP_,  L0_,  Ll_,  NL_]  :=  Module [ 

{XO,  XI,  nx,  YO,  Yl,  ny,  Ftix,  Mhue,  Clevels,  XF}, 

XO  =  XP [ [ 1 ] ] ;  XI  =  XP [ [ 2 ] ] ;  nx  =  XP [ [ 3 ] ] ; 

YO  =  XP [ [ 4 ] ] ;  Yl  =  XP [ [ 5 ] ] ;  ny  =  XP [ [ 6 ] ] ; 

PP  =  XP [ [ 7 ] ] ; 

Mhue[h_]  =  Hue[l,  0,  1] ; 

Ftix  =  {{{1,  XO],  {0.75  1  +  0.25  nx,  0.75  XO  +  0.25  XI}, 
{0.5  1  +  0.5  nx,  0.5  X0  +  0.5  XI}, 

{0.25  1  +  0.75  nx,  0.25  X0  +  0.75  XI},  {nx,  XI}}, 
{{1,  Y0},  {0.5  1  +  0.5  ny,  0.5  Y0  +  0.5  Yl}, 

{ny,  Yl}},  None,  None}; 

Clevels  =  Log[Lseq[L0,  LI,  NL]]; 

Print [N [Round [1  Exp [Clevels] ]  /  1] ] ; 

XF  =  Lis tContourPlot [Transpose [Log [PP] ] , 

AspectRatio  ->  1/2,  ColorFunction  ->  Mhue, 
FrameTicks  ->  Ftix,  Contours  ->  Clevels, 

Contour Smoothing  ->  32]; 

XF] 


A.3  Example  Use.  The  graphics  in  this  report  were  produced  by  the  following 
commands.  First,  create  the  numerical  arrays.  With  nn=  200,  the  functions  are  evaluated 


on  a  200  x  100  grid.  This  takes  a  while, 
lower  resolution  results. 

nn  =  200; 

XP1  =  Ptable [20,  nn]  ; 

XP2  =  Ptable[10,  nn] ; 

XP3  =  Ptable [5,  nn]  ; 

XP4  =  Ptable [2,  nn] ; 

XP5  =  Ptable [1,  nn] ; 

Then,  construct  the  graphs. 

XFl  =  Pshow [XP1,  2,  20,  7] 

XF2  =  Pshow [XP2 ,  5,  100,  7] 

XF3  =  Pshow [XP3,  10,  500,  7] 

XF4  =  Pshow [XP4,  20,  1500,  7] 

XF5  =  Pshow [XP 5 ,  50,  3000,  7] 


Smaller  values  of  nn  can  be  used  for  quicker. 
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Finally,  export  the  graphics  files. 


IS  =  {600,  600}; 
Display [ "xfl .ps" ,  XFl, 
Display [ "xf 2 .ps" ,  XF2 , 
Display ["xf 3 .ps" ,  XF3 , 
Display [ "xf 4 .ps" ,  XF4 , 
Display [ "xf 5 .ps" ,  XF5, 


EPS",  ImageSize  ->  IS] ; 
EPS",  ImageSize  ->  IS]; 
EPS",  ImageSize  ->  IS]; 
EPS",  ImageSize  ->  IS]; 
EPS",  ImageSize  ->  IS]; 
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6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to  organization, 
technical  content,  format,  etc.) _ _ _ • _ 


Organization 

CURRENT  Name  E-mail  Name 

ADDRESS  _ _ _ 

Street  or  P.O.  Box  No. 

City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the  Old 
or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 

(DO  NOT  STAPLE) 


